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Abstract 

A new method is presented to determine the gluon density in the proton from jet pro- 
duction in deeply inelastic scattering. By using the technique of Mellin transforms not only 
for the solution of the scale evolution equation of the parton densities but also for the eval- 
uation of scattering cross sections, the gluon density can be extracted in next-to-leading 
order QCD. The method described in this paper is, however, more general, and can be 
used in situations where a repeated fast numerical evaluation of scattering cross sections 
for varying parton distribution functions is required. 
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1 Introduction 



One of the main goals of experiments at the electron-proton collider HERA is the precise deter- 
mination of the gluon density f g / p (^,fi 2 ) in the proton for various gluon momentum fractions £ 
and factorization scales \i. In addition to the indirect method of extracting f g / p from the scaling 
violation of the structure function F2, direct methods such as heavy quark and jet production 
have been studied. 

In the QCD-improved parton model, the electron-proton scattering cross section a is gener- 
ically given by a convolution of process-independent parton densities f q / p for (anti-) quarks and 
fg/p for gluons with corresponding mass-factorized parton- level cross sections a q and a g : 



= J d ^ [fv/p^'^Vqitv 2 ) + fg/p(£,t* 2 )<7g(€,fi l 



Besides the indicated dependence on £ and fx, a q and a g also depend on other variables such 
as the absolute electron four-momentum transfer squared Q 2 and the momenta of the outgoing 
partonsQ. 




Figure 1: Generic Feynman diagrams for the leading-order processes of QCD Compton scattering 
(a) and photon- gluon fusion (b), and an example for a diagram corresponding to a next-to-leading 
order real correction (c). 

In leading order (LO), the prescription for the extraction of f g / p from jet cross sections in 
deeply inelastic scattering reactions is very intuitive, because a q can be identified with the parton- 
model cross section o~c for the so-called QCD Compton scattering process (fig. |l]a), and o~ g with the 
parton-model cross section o~p for the photon- gluon fusion reaction (fig. |I]b). Explicit expressions 
can be found in [0J. Experimentally the outgoing partons from the hard scattering reactions are 
identified with jets. The QCD Compton scattering and photon-gluon fusion reactions lead to 
(2+l)-jet final states, where the notation accounts for the two outgoing jets from the hard 
scattering process and the jet in the proton fragmentation region. The calculated contribution 

1 Here and in the following we do not explicitly display the dependence on the renormalization scale. 



<7 c%+i = I d£ fq/piC, H 2 ) a c(0 from Compton scattering can be subtracted from the measured 
cross section 



and thus f g / P {^,f^ 2 ) can be determined in LO by a direct unfolding, since in this case £ can 
be expressed in terms of measurable quantities as £ = xb (1 + $/Q 2 ), where xb is the Bjorken 
scaling variable and s is the invariant mass squared of the system of the two current jets. An 
analysis based on this principle has recently been presented by the HI Collaboration ||. 

In next-to-leading order (NLO) this simple picture is destroyed. Aside from the virtual 
corrections to the Born processes in figs. |T]a and b, real corrections have to be added; diagrams 
of the type shown in fig. |l]c can also lead to (2+l)-jet configurations: if the gluon g attached 
to the outgoing quark is soft or collinear to the quark, the diagram constitutes a correction to 
the photon-gluon fusion process. If, on the other hand, this gluon is hard and the outgoing 
antiquark q is soft or collinear to the incoming gluon g', then this configuration can be said to 
be a correction to the QCD Compton scattering reaction. In the latter case, the collinear or soft 
antiquark forms a jet with the proton remnant r, and the cross section has to be integrated over 
all momenta of the antiquark according to a specific jet definition scheme. Collinear singularities 
that do not cancel against corresponding singularities from the virtual corrections have to be 
absorbed into renormalized parton densities. Depending on the factorization scheme chosen, 
finite subtracted pieces remain. The factorization theorems of perturbative QCD guarantee that 
the cross section can be written in the form of eq. ([!]). However, beyond the leading order, 
the arbitrary momentum of collinear partons renders the variable £ unobservable, because the 
mass-factorized parton-level cross sections are in general distributions, not regular functions, 
and the simple and straightforward method described above can therefore not be applied. A 
physical consequence is that the distinction between the QCD Compton and photon-gluon fusion 
processes becomes meaningless. Related to this is the fact that quark and gluon densities mix in 
the Altarelli-Parisi scale evolution. 

A determination of the gluon density in NLO is very desirable. In LO the partonic cross 
sections a q and a g (in short denoted by cr^) do not depend on /z, as already indicated in eq. (H), 
and the fi/ p (as short-hand for f q / p and f g / p ) are the solutions of the LO Altarelli-Parisi evolution 
equation, where the leading logarithmic terms in the scale /i are summed up. In any finite order 
of perturbation theory, the scattering cross section a depends explicitly on the factorization scale 
/i, this scale dependence being due to uncalculated higher-order terms. The scale dependence is 
particularly strong in the LO case, because there no compensation can take place between f^/ p 
and Oi. To a great extent this problem is, for many processes, reduced in NLO, where explicit 
terms ~ ln/i 2 in Oi compensate the /i-dependence of fi/ p such that the variation is of higher 
order in the strong coupling constant a s . For reliable theoretical predictions, a NLO analysis of 
scale-dependent quantities is therefore mandatory. 

The only way to achieve a direct NLO determination of f g / p is to parametrize the function 




(2) 



fg/ p at a given scale /i , to evolve it then to a value of \i where the cross section is measured, say 
/i = Q, and to fit the parameters of f g / p with respect to suitable infrared safe observables, e.g. 
the (2+l)-jet cross section in various bins of xb- A severe practical problem is that the cross 
section a has to be evaluated repeatedly for every choice of parameters for f g / p . Monte Carlo 
methods allow the application of arbitrary cuts on final-state particle momenta, as is necessary in 
order to take detector acceptance cuts properly into account, but these methods are prohibitively 
slow. A fast numerical method for the repeated application of this procedure is indispensable, 
and will be developed in this paper. 

The paper is organized as follows. The new method is formally derived in Section 0. Details 
of the Mellin transform relevant to the application of the method are discussed in Section |3|. 
Finally, an explicit numerical example is given in Section ^ for the case of jet cross sections in 
deeply inelastic electron-proton scattering, where it is shown that the method is operational in 
practice when realistic acceptance cuts are taken into account. The paper closes with a short 
summary. 



2 The Mellin Transform Technique 
for Non- Factorizing Cross Sections 

The Mellin transform technique allows for a quick numerical evaluation of integrals of the form 

(3) 



S(ar B ) = y/i/ P (O ff i ^ x b^ 



in the case where is independent of its second argument xb, on the basis of the moments 
defined by 

F n = f 1 — x n F(x) (4) 
Jo X 

for an arbitrary function F and (complex) n. The moments of the function E are then given by 

^-"n fi/p,n &i,n- (5) 

The functional dependence of E can be recovered from the moments E n by an inverse Mellin 
transform. An expression of the form of eq. (|f) will be called to be of the factorizable type if the 
only dependence on xb in the arguments of <jj is via xb/C,- In the application which we have 
in mind, fi/ p is a parton density, whereas <Tj is an expression for a mass-factorized parton- level 
scattering cross section. In general, acceptance cuts and non-factorizable jet algorithms (cf. the 
discussion in |3|, |J) introduce an explicit dependence of <jj on xb- Moreover, the expression for 
T,(xb) is integrated over a certain range of xb- This might suggest that the Mellin transform 
technique cannot be applied. However, this is not the case. In the following we outline a method 
that allows the use of this technique. 



The cross section differential in xb can be written in the form 



= f — fi/ p {€, /x 2 ) Ui 



Jx B t, 



( 




(6) 



where 




) 




dT &i — ,x B ,T,/i 2 




(7) 



Here we have made the dependence on the factorization scale /i explicit. The set T of variables 
contains all other integration variables besides xb, i-e. the other electron variables including Q 2 
and the momenta of the outgoing partons (or jets); and V XB is the phase space region over which 
these variables are integrated. Thus, V XB includes all acceptance and jet cuts, as well as the range 
in Q 2 for the specified xb- <5% is the cross section differential in all variables including T, whereas 
a - , is the integration kernel to be convoluted with jq v to yield S. The explicit dependence on xb 
of <Ji(xB/£,XB,T, /i 2 ) and a,i(xB / '£,xb, IJ 2 ) (not just via the ratio Xb/C) is a consequence of the 
(finite) factorization breaking terms^. 

It is assumed that the factorization scale /i is chosen independently from the variables T; 
in particular, // should not depend on the integration variable Q 2 , which is integrated over a 
certain range [Qo,<3 2 ]. // it may be set to some intermediate value. In any case, this is not a 
strong restriction, because the parton densities depend only logarithmically on the factorization 
scale, and moreover the scale dependence is compensated by a corresponding term in the mass- 
factorized parton- level scattering cross section, so that the change is of higher order in a s . For 
simplicity of notation, we drop the argument /i 2 in the following. 

Now let (for a fixed Q 2 -bin) ai, . . . , a& be the experimental boundaries of the intervals in the 
variable xb for which the cross sections are measured. To proceed, we define 



In a fit of the function fi/ p , the integrals in eqs. ©, ([?]), (§) have to be evaluated repeatedly. In 
general, the xb, £ and T-integrations are performed by a time-consuming Monte Carlo integration 
in order to implement all cuts. In particular, the (xb, Q 2 )-plane is divided into several bins which 
do not change during the fitting procedure. The method requires that a certain set of moments 
is calculated for every bin in the (xb, Q 2 )-plane. The steps to be followed to determine the S a 
are: 

2 It should be kept in mind that the infinite terms related to infrared and collinear singularities still factorize in 
the standard form so that universal parton densities can be defined. This property is not spoiled by acceptance 
cuts and non- factorizing jet definition schemes. 





The integral over a specified interval [ai, in xb is then simply given by 




(9) 



i 



Define a function 



{ J ( X B 

/ dx B Oi — - 

< J a \a/u 



ib , if u > a 



h a {u) = \ Ja \a/u /' (10) 

if u < a 



and its moments in the variable w 



= / — U n h a {u). (11) 
JO M 



It is easy to prove that an explicit expression for h an is 

Kn = r ^ l f (?) ^(f'^)- ^ 

The /i an are thus the S a with the parton density fi/ p (£) replaced by (a/£) n . They can be 
determined numerically by means of a Monte Carlo integration. In general, for complex n, 
the quantity (a/£) n has to be split into its real and imaginary part. 

Define 

^ = [ f /</ p (fl h (£j (13) 
and determine the moments of E a b with respect to the variable a: 

S n6 = / — a n E ab . (14) 
Jo a 

Obviously, S a = S aa . The key relation of our method is 

^nfe = fi/p,n hbni (15) 

and can be proved in the following way: 



/:fw>»(i) 

1 df dw 
o £ Jo u 

— fi/p,n hbn- (16) 

For a given parametrization of /j/ p in terms of its moments fi/ pn the cross section S a can 
be determined by forming the moments T, na = fi/ p>n h an and a subsequent inverse Mellin 
transform in the variable n, evaluated at a. 



The Mellin transform method in the case of non-factorizing cross sections introduces the 
inconvenience that the moments h an have to be determined for every interval boundary a sep- 
arately. Due to the large number of repeated cross section evaluations in the fitting procedure, 
however, this method is far more efficient than a direct integration of the integrals in eqs. @, 
® f° r every parametrized parton density. 



3 From Parton Moments to Observables 



Let us now consider the inverse transformation of the moments given by eq. (f|), which is a special 
case of the general Mellin transformation for functions F(x) vanishing identically at x > 1. If 
F(x) is piecewise smooth for x > 0, the corresponding Mellin inversion reads 

I rc+ioo 

F{x) = — : / dnx~ n F n , (17) 



2m 



C—lOO 



where the real number c has to be chosen such that Jq 1 dxx c ~ 1 F(x) is absolutely convergent 
IJ. Hence c has to lie to the right of the rightmost singularity n max of F n . The contour of the 
integration in eq. (|17D is displayed in fig. § and denoted by C . Also shown is a deformed route Ci, 
yielding the same result as long as no singularities rii of F n are enclosed by Cq — Ci. For example, 
for the LO and NLO evolution of structure functions, the n» are real with m < n max < c, and 
this requirement is fulfilled automatically. 
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Figure 2: Integration contours of the Mellin inversion in eq. (\F/[), leading to the inversion formu- 
lae of eqs. filHl ) and fcU\ ) for the routes Co and C\, respectively. The crosses schematically denote 
the singularities of F n . 



It is useful to rewrite eq. (|TT| ) as an integration over a real variable. We are concerned with 
functions obeying F* = F n *, where '*' denotes the complex conjugation. Then it is easy to show 
that eq. (|T7|) yields, for the contour characterized by the abscissa c and the angle in fig. |2|: 

1 



F(x) 



7T JO 



dz Im 



exp (i0)aT c - zexpW) F, 



n=c+z exp (iq 



It is obvious from the discussion given above that the integral does not depend on c and 0. 
However, for an efficient numerical evaluation a suitable choice of these parameters is very useful. 
For example, it is advantageous to choose <fi > tt/2 in case F n is a known analytical function, 



especially if this function does not fall off very rapidly for \n\ — > oo. The factor exp (z log ~ cos <p\ 
then introduces an exponential dampening of the integrand (which rapidly oscillates at small x) 
with increasing z, thereby allowing for a smaller upper limit z max in the numerical implementation 
of eq. (|18|). This procedure has been employed for the inversion of moments of parton densities 
and structure functions for the proton and the photon, e.g. in || [7[] and [§], respectively. 

In general, however, the moments of the partonic cross section can only be calculated numer- 
ically using eq. because no analytic continuation to small Ren, where the integral does not 
exist, is at our disposal. Likewise, in our case these moments are given by eq. (|12D and do not 
behave uniformly for |n| — ► oo. Especially, they grow exponentially along C\. Therefore, we will 
use the 'textbook contour' Co in the following and, with = n/2, eq. (|1^) simplifies to 



Fix) = - dzRe 

7T JO 



oo 



n=c+iz 



(19) 



We have applications in mind where F n = fi/ Ptn h an , see eq. flT5p, and the numerical evaluation 
of the moments h an in eq. (|12|) is very time-consuming. Taking a different upper limit z max of 
the numerical z-integration or number of points for the integral at each step in the integration 
process is practically unfeasible in such a case. Instead, we want to fix z max at a value as small 
as possible in order to allow for an evaluation of eq. (|19|) with a rather small number of fixed 
moments. 

In this context, it is useful to consider that part of the n-dependent inversion integrand, of 
which the analytical continuations are known together with the behaviour at large n, namely the 
parton densities and their evolution. Inspection of eq. (|T2"| ) suggests that h an does not rise strongly 
with z along C - Hence the large-n behaviour of the parton densities /i/ Pi „ can be employed to 
estimate the convergence of the complete integral in eq. (|19"D with F n given by eq. fli~5|) . A typical 
ansatz for the parton distribution functions of the proton at some reference scale /iQ, denoted by 

/»/p(0) is g iven b y M> S 



e/i/p(0=^r(l-fl'(l + 7)/e + --O • (2°) 

The coefficients j3 can be estimated roughly by their counting rule values, e.g. /3 va i ~ 3, /3 g \ ne ~ 5. 
The Mellin transform of eq. (|20| ) reads simply 



f i/p>n = A [B{a + n - 1, /3 + 1) + jB{a + n- 1/2, (5 + 1) + . . .] (21) 
with the Euler Beta-function B. If (5 is a positive integer, then this equation simplifies to 

fi/P,n = r— Tw , A f' ~ —„ 7\ + • • • = 0(l/n? +1 ) for n > oo . (22) 

[a + n — 1) (a + n) . . . [a + n + p — 1) 

The evolution of these input moments is known analytically for arbitrary complex n 0], and 
the kernel Kij^fi 2 , /Zq) in 

h/pAv 2 ) = KijAv 2 ' iA) fj/pAvl) ( 23 ) 



generally leads to a slightly faster decrease of fi/ p>n for large n at /i 2 > /!q. Hence a fall-off like 
l/n 4 can be safely used to estimate the practically required upper limit in eq. flT9"|). 

For this purpose, we have numerically determined the upper limit z max sufficient to reach a 
1% accuracy of the Mellin inversion for the toy function F(x) = — x) 3 . The results are 

displayed in table |I| for selected values of x. The rightmost pole is at n max = 1, and we have 
chosen c = 1.5. The larger c is, the more F n=c+iz is flattened. Hence too large a value of c leads 
to an undesired rise of z max at small x, where F n is integrated after multiplication with a rapidly 
oscillating function in eq. ([19]). Practically, c — n max ~ 0.5 works well for all x- values of interest 
here, implying c ~ 1.8 for realistic small-^ parton densities |7], |J. 



X 


10" 4 


10" 3 


0.01 


0.03 


0.1 


0.3 


0.8 


Zmax 


6.0 


5.5 


4.5 


4.0 


3.5 


5.0 


8.0 



Table 1: Upper limits z max numerically sufficient for a 1% accuracy of the Mellin inversion from 
eq. ( \TQ ) for the function F(x) = — x) 3 with c = 1.5. 



The integral eq. (|T9|), truncated at z max , can now be performed by using a sufficiently large 
number of fixed support points, e.g. by a sum of 8-point Gaussian quadratures, see fTTJI for 
the weights and support points. In this way, everything except for the input-5-functions varied 
in a fit of the parton distribution functions is fixed; especially the time-consuming part of the 
kernels in eq. (|23|) and the moments h an of the partonic cross sections from eq. (|i~2"D can be 
determined once and then used unchanged in the calculation of physical observables for various 
parton densities. 



4 Application to Jet Physics at HERA 

To illustrate how the Mellin transform method can be used to fit the gluon density f g / p , the gluon- 
induced (2+l)-jet cross sections were calculated in several bins for HERA energies of 820 GeV 
protons and 27.6 GeV electrons. Quark contributions were set to zero explicitly in the parton 
distribution function to reduce the number of moments needed for this case study. 



The program PROJET Jl2[ based on the NLO matrix elements from |T3| , [T4]| was used, this 
allows to calculate jet cross sections in LO and NLO in the modified JADE scheme defined in 
the following way [151]: 



• Define a precluster of longitudinal momentum p r given by the missing longitudinal momen- 
tum of the event. 

o 



• Apply the JADE cluster algorithm [|HJ to the set of momenta {pi, . . . ,p n ,p r }, where 
Pi, . . . ,p n are the momenta of the hadrons visible in the detector. The resolution crite- 
rion is Sij = IpiPj > y cu tM 2 . Here M 2 is a mass scale and y cut is the resolution parameter. 



In the case of a theoretical calculation, p r is directly given by the momentum fraction of the 
proton not carried by the incident parton, and pi, . . . ,p n are the momenta of the partons in the 
final state. In the following, we choose W 2 , the squared total hadronic energy, as the mass scale 
M 2 , since the proton remnant is included in the jet definition. 



The integration routine used in PRO JET is VEGAS |17|, |l8 |. As is desirable for an exper- 
imental measurement, the phase space was binned in Q 2 and xb according to eq. (||); the bins 
are given in tables ^ and |3|. In addition, the following typical HI detector cuts were applied, for 



which the motivation is explained in detail in [19|: 



The invariant mass squared of the hadronic system W 2 was required to be larger than 
5000 GeV 2 . 

The jet resolution cut y cut was set to 0.02. Lowering this value significantly below 0.01 
causes NLO corrections to dominate and leads to unphysical cross sections. It is important 
to note that £ > y cut as a consequence of the applied modified JADE algorithm. The region 
£ > 0.01 is however very interesting [2(J for a precise determination of f g / p , see also |2l| . 



The jets were required to lie in the polar angle range of 10° < 6j et < 145°. 

For bins with Q 2 < 100 GeV 2 , the scattered electron had to have an energy of Ep > 14 GeV 
and the polar angle had to lie within the range of 160° < Qy < 172.5°. 

In the bins with Q 2 > 100 GeV 2 , the scaled photon energy y in the proton rest system had 
to be y < 0.7 and the scattered electron was required to have 10° < Qy < 148°. 



In this list, angles and energies are defined in the laboratory frame, and angles are given with 
respect to the direction of the incoming proton. For each bin, 32 complex Mellin moments were 
calculated according to the prescription described in Section 0, cf. eq. (|T2"D. In all calculations, 
a s was computed to second order, and the NLO gluon distribution function of [[/]] was employed. 

A good convergence of the numerical calculations was found for c = 1.8, <fi = ir/2 and 
Zmax = 9, with a higher density of support points at lower z, as the influence is greatest there. 
For comparison, the cross section was also calculated directly, see eqs. (||), (0), (H). After inverting 
the product of the hard subprocess and evolved gluon density moments at the average Q 2 , the 
results were found to coincide at the per cent level. The detailed results can be found in tables || 
and |||. In most bins, convergence was reached at z max = 3 (corresponding to 16 moments), the 
additional moments were used for safety. The convergence of the LO cross section was much faster 
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Q 2 [GeV 2 ] 


1 A — 4 

10 


. . . 1 


i n— 3 
10 


. . . 1 


i n— 2 
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. . . 1 
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10 
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14 


62.80 
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28.09 
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14 .. 


18 


70.64 


69.72 


50.95 


49.97 
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25 


85.82 
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40 


109.9 
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101.9 


101.1 










40 .. 


100 






123.8 


124.4 


14.51 


14.43 






100 .. 


300 






31.96 


32.18 


14.69 


14.76 






300 .. 


700 






28.97 


29.23 


25.18 


25.42 






700 .. 


. 4000 










10.22 


10.12 


0.96 


0.93 



Table 2: Comparison of cross sections with LO matrix elements 3 (in [pb]) obtained by integrating 
directly (left columns) or using the Mellin transform method (right columns). 





XB 


Q 2 [GeV 2 ] 


10~ 4 


...1 


10~ 3 
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58.48 


57.25 


26.60 


26.00 
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18 


66.57 


65.90 


47.22 


46.69 
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82.48 


81.65 


67.99 
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100.4 
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125.6 


14.07 
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300 






34.86 


34.52 


15.51 


15.31 






300 .. 


700 






31.34 


31.51 


27.01 


27.19 






700 .. 


. 4000 










11.18 


11.19 


0.99 


0.97 



Table 3: Comparison of NLO cross sections (in [pb]) obtained by integrating directly (left columns) 
or using the Mellin transform method (right columns). 



than in the NLO for a given number of support points in VEGAS, the LO integration 

is more accurate due to the simpler integration kernel. The method works well for both LO and 
NLO. 

The number of points in the Monte Carlo integration was chosen such that the error returned 



by VEGAS was less than 1%. This number is, however, only a rough estimate |L7|> [18|| , and 
the achieved accuracy was studied by repeating the calculation for different random number 
generator seeds. The direct integrations performed here had a statistical variation of 2-3%. The 
partonic cross section from the Mellin transform method is implicitly integrated repeatedly by 
the calculation of the moments, which smoothes out statistical variations. The results were found 



3 Here, 'LO' means that the matrix elements were calculated in LO, but a s and the parton distribution functions 
in NLO to facilitate a comparison with the results of table ||. For a physically meaningful comparison of the LO 
with the NLO, a s and the parton distribution functions should be calculated in LO, if they are used in conjunction 
with the LO matrix elements. 



i n 



to be more stable than the direct integration, which varied around the result obtained by the 
moment inversion. Even drastic errors of single moments or setting single moments to zero could 
be tolerated and led to a reproducible result. We conclude that this method is numerically very 
stable and that the accuracy is of the order of 1%. Increasing the accuracy requires increasing the 
number of support points for the integration, which would result in a dramatic increase in CPU 
timeQ. One has to keep in mind that an additional error source arises from the Mellin transform 
method, as for each experimental bin in x one has to calculate the difference of the cross sections 
depending on the bin boundaries in eq. (Q), leading to error propagation. A strategy for bin 
optimization is under study. 

5 Summary 

We have outlined a new method, based on the Mellin transform for the fast evaluation of the 
convolution of a parton-level cross section with a parton density, which works for cross sections 
not showing a simple factorization behaviour. The method can be the basis for a fit of the gluon 
density from experimental data for the (2+l)-jet cross section in deeply inelastic electron-proton 
scattering, but it is also suitable for more complicated observables. It has been explicitly shown 
by a numerical study using realistic experimental cuts, that the method works in practice with 
a sufficiently high accuracy. 
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1 Introduction 



One of the main goals of experiments at the electron-proton collider HERA is the precise deter- 
mination of the gluon density / g / p (£, fi 2 ) in the proton for various gluon momentum fractions £ 
and factorization scales fi. In addition to the indirect method of extracting f g / p from the scaling 
violation of the structure function F 2 , direct methods such as heavy quark and jet production 
have been studied. 

In the QCD-improved parton model, the electron-proton scattering cross section a is gener- 
ically given by a convolution of process-independent parton densities f q / p for (anti-)quarks and 
f g / p for gluons with corresponding mass-factorized parton-level cross sections a q and a g : 

° = j^[f^^ 2 )^^ 2 ) + f 9 / P ^^>M^ 2 )]- (1) 

Besides the indicated dependence on £ and /i, a q and a g also depend on other variables such 
as the absolute electron four-momentum transfer squared Q 2 and the momenta of the outgoing 
partons 1 . 




Figure 1: Generic Feynman diagrams for the leading-order processes of QCD Compton scattering 
(a) and photon- gluon fusion (b), and an example for a diagram corresponding to a next-to-leading 
order real correction (c). 

In leading order (LO), the prescription for the extraction of f g / p from jet cross sections in 
deeply inelastic scattering reactions is very intuitive, because a q can be identified with the parton- 
model cross section o~c for the so-called QCD Compton scattering process (fig. la), and a g with the 
parton-model cross section ap for the photon-gluon fusion reaction (fig. lb). Explicit expressions 
can be found in [1]. Experimentally the outgoing partons from the hard scattering reactions are 
identified with jets. The QCD Compton scattering and photon-gluon fusion reactions lead to 
(2+l)-jet final states, where the notation accounts for the two outgoing jets from the hard 
scattering process and the jet in the proton fragmentation region. The calculated contribution 

1 Here and in the following we do not explicitly display the dependence on the renormalization scale. 
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a c°2+i = I fq/p{ii Z^ 2 ) a c{C) from Compton scattering can be subtracted from the measured 
cross section 

^2° = f*t [/ s/ p(e,M 2 )Mo + f,/p(t,p)Mtj\ , ( 2 ) 

and thus / fl / p (£, fi 2 ) can be determined in LO by a direct unfolding, since in this case £ can 
be expressed in terms of measurable quantities as £ = (1 + s/Q 2 ), where is the Bjorken 
scaling variable and s is the invariant mass squared of the system of the two current jets. An 
analysis based on this principle has recently been presented by the HI Collaboration [2]. 

In next-to-leading order (NLO) this simple picture is destroyed. Aside from the virtual 
corrections to the Born processes in figs, la and b, real corrections have to be added; diagrams 
of the type shown in fig. lc can also lead to (2+l)-jet configurations: if the gluon g attached 
to the outgoing quark is soft or collinear to the quark, the diagram constitutes a correction to 
the photon-gluon fusion process. If, on the other hand, this gluon is hard and the outgoing 
antiquark q is soft or collinear to the incoming gluon g' ', then this configuration can be said to 
be a correction to the QCD Compton scattering reaction. In the latter case, the collinear or soft 
antiquark forms a jet with the proton remnant r, and the cross section has to be integrated over 
all momenta of the antiquark according to a specific jet definition scheme. Collinear singularities 
that do not cancel against corresponding singularities from the virtual corrections have to be 
absorbed into renormalized parton densities. Depending on the factorization scheme chosen, 
finite subtracted pieces remain. The factorization theorems of perturbative QCD guarantee that 
the cross section can be written in the form of eq. (1). However, beyond the leading order, 
the arbitrary momentum of collinear partons renders the variable £ unobservable, because the 
mass-factorized parton-level cross sections are in general distributions, not regular functions, 
and the simple and straightforward method described above can therefore not be applied. A 
physical consequence is that the distinction between the QCD Compton and photon-gluon fusion 
processes becomes meaningless. Related to this is the fact that quark and gluon densities mix in 
the Altarelli-Parisi scale evolution. 

A determination of the gluon density in NLO is very desirable. In LO the partonic cross 
sections a q and a g (in short denoted by <t;) do not depend on /i, as already indicated in eq. (2), 
and the fi/ p (as short-hand for f q / p and f g / p ) are the solutions of the LO Altarelli-Parisi evolution 
equation, where the leading logarithmic terms in the scale fi are summed up. In any finite order 
of perturbation theory, the scattering cross section a depends explicitly on the factorization scale 
/i, this scale dependence being due to uncalculated higher-order terms. The scale dependence is 
particularly strong in the LO case, because there no compensation can take place between fi/ p 
and To a great extent this problem is, for many processes, reduced in NLO, where explicit 
terms ~ In fi 2 in <t; compensate the //-dependence of fi/ p such that the variation is of higher 
order in the strong coupling constant a s . For reliable theoretical predictions, a NLO analysis of 
scale-dependent quantities is therefore mandatory. 

The only way to achieve a direct NLO determination of f g / p is to parametrize the function 
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fg/p at a given scale /i , to evolve it then to a value of fi where the cross section is measured, say 
fi = Q, and to fit the parameters of f g / p with respect to suitable infrared safe observables, e.g. 
the (2+l)-jet cross section in various bins of xb- A severe practical problem is that the cross 
section a has to be evaluated repeatedly for every choice of parameters for f g / p - Monte Carlo 
methods allow the application of arbitrary cuts on final-state particle momenta, as is necessary in 
order to take detector acceptance cuts properly into account, but these methods are prohibitively 
slow. A fast numerical method for the repeated application of this procedure is indispensable, 
and will be developed in this paper. 

The paper is organized as follows. The new method is formally derived in Section 2. Details 
of the Mellin transform relevant to the application of the method are discussed in Section 3. 
Finally, an explicit numerical example is given in Section 4 for the case of jet cross sections in 
deeply inelastic electron-proton scattering, where it is shown that the method is operational in 
practice when realistic acceptance cuts are taken into account. The paper closes with a short 
summary. 

2 The Mellin Transform Technique 
for Non-Factorizing Cross Sections 

The Mellin transform technique allows for a quick numerical evaluation of integrals of the form 

S(aj B ) = J y/i/plO^lyi^j (3) 

in the case where <t; is independent of its second argument xb-, on the basis of the moments 
defined by 

F n = C — x n F(x) (4) 

JO X 

for an arbitrary function F and (complex) n. The moments of the function S are then given by 

fi/p,n &i,n- (5) 

The functional dependence of S can be recovered from the moments S n by an inverse Mellin 
transform. An expression of the form of eq. (3) will be called to be of the factorizable type if the 
only dependence on xb in the arguments of <t; is via xbH- In the application which we have 
in mind, is a parton density, whereas <t; is an expression for a mass-factorized parton-level 
scattering cross section. In general, acceptance cuts and non-factorizable jet algorithms (cf. the 
discussion in [3, 4]) introduce an explicit dependence of <t; on xb- Moreover, the expression for 
Ti(xb) is integrated over a certain range of xb- This might suggest that the Mellin transform 
technique cannot be applied. However, this is not the case. In the following we outline a method 
that allows the use of this technique. 
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The cross section differential in xb can be written in the form 



where 




Here we have made the dependence on the factorization scale fi explicit. The set T of variables 
contains all other integration variables besides xb, i-e. the other electron variables including Q 2 
and the momenta of the outgoing partons (or jets); and V XB is the phase space region over which 
these variables are integrated. Thus, V XB includes all acceptance and jet cuts, as well as the range 
in Q 2 for the specified xb- &i is the cross section differential in all variables including T, whereas 
(Ji is the integration kernel to be convoluted with fi/ p to yield S. The explicit dependence on xb 
of p^xbUi XB, fi 2 ) and cr^xs/i, xb, fi 2 ) (not just via the ratio xbH) is a consequence of the 
(finite) factorization breaking terms 2 . 

It is assumed that the factorization scale fi is chosen independently from the variables T; 
in particular, fi should not depend on the integration variable Q 2 , which is integrated over a 
certain range [£?q,£?i]. f 1 ^ mav be set to some intermediate value. In any case, this is not a 
strong restriction, because the parton densities depend only logarithmically on the factorization 
scale, and moreover the scale dependence is compensated by a corresponding term in the mass- 
factorized parton-level scattering cross section, so that the change is of higher order in a s . For 
simplicity of notation, we drop the argument fi 2 in the following. 

Now let (for a fixed Q 2 -b'm) a\, . . . , a& be the experimental boundaries of the intervals in the 
variable xb for which the cross sections are measured. To proceed, we define 

T, a = f dx B T,(xb). (8) 

J a 

The integral over a specified interval [a;,a; + i] in xb is then simply given by 

/ dx B T,(x B ) = S ai - S ai+1 . (9) 

J Cti 

In a fit of the function the integrals in eqs. (6), (7), (8) have to be evaluated repeatedly. In 
general, the xb-, £ and T-integrations are performed by a time-consuming Monte Carlo integration 
in order to implement all cuts. In particular, the (xb-, <5 2 )-plane is divided into several bins which 
do not change during the fitting procedure. The method requires that a certain set of moments 
is calculated for every bin in the {xb-, <5 2 )-plane. The steps to be followed to determine the S a 
are: 

2 It should be kept in mind that the infinite terms related to infrared and collinear singularities still factorize in 
the standard form so that universal parton densities can be defined. This property is not spoiled by acceptance 
cuts and non-factorizing jet definition schemes. 
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• Define a function 



r /u i / x b 

/ dx B (Ti — — 

Ja \a/u 



xb , if u > a 



h a (u) = { Ja \a/u J (10) 

0, if ti<d 



and its moments in the variable u 

du 



r du 

h an = / — u n h a {u). (11) 
Jo u 

It is easy to prove that an explicit expression for h an is 

han = la dXB Cj ^(T'^)- (12) 

The h an are thus the S a with the parton density fi/ p (i) replaced by (a/£) n . They can be 
determined numerically by means of a Monte Carlo integration. In general, for complex n, 
the quantity (a/^) n has to be split into its real and imaginary part. 



• Define 

Sa6 = £ f MOhb {f) (13) 
and determine the moments of S a & with respect to the variable a: 

t nb = f 1 —a n H ab . (14) 
Obviously, S a = S aa . The key relation of our method is 

^nb = fi/p,n,hbni (^) 

and can be proved in the following way: 

^ = /:^-/:f^«^(?) 

= t %■ r mo t — « n h(u) 

Jo c, Jo u 

= fi/p,nhbn- (16) 

• For a given parametrization of in terms of its moments fi/ p>n the cross section S a can 
be determined by forming the moments S na = fi/ Pin h an and a subsequent inverse Mellin 
transform in the variable n, evaluated at a. 



The Mellin transform method in the case of non-factorizing cross sections introduces the 
inconvenience that the moments h an have to be determined for every interval boundary a sep- 
arately. Due to the large number of repeated cross section evaluations in the fitting procedure, 
however, this method is far more efficient than a direct integration of the integrals in eqs. (6), 
(7), (8) for every parametrized parton density. 
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3 From Parton Moments to Observables 



Let us now consider the inverse transformation of the moments given by eq. (4), which is a special 
case of the general Mellin transformation for functions F(x) vanishing identically at x > 1. If 
F(x) is piecewise smooth for x > 0, the corresponding Mellin inversion reads 

1 /»c+ioo 

F(x) = — dnx- n F n , (17) 

Z7TZ J c—ioo 

where the real number c has to be chosen such that Jq dx x c ~ 1 F(x) is absolutely convergent 
[5]. Hence c has to lie to the right of the rightmost singularity n max of F n . The contour of the 
integration in eq. (17) is displayed in fig. 2 and denoted by Co- Also shown is a deformed route C\, 
yielding the same result as long as no singularities rt{ of F n are enclosed by Cq — C\. For example, 
for the LO and NLO evolution of structure functions, the rt{ are real with rt{ < n max < c, and 
this requirement is fulfilled automatically. 



\ 

\ 
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X X X ) 
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V / 
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c Ren 



Figure 2: Integration contours of the Mellin inversion in eq. (17), leading to the inversion formu- 
lae of eqs. (18) and (19) for the routes Co and C\, respectively. The crosses schematically denote 
the singularities of F n . 

It is useful to rewrite eq. (17) as an integration over a real variable. We are concerned with 
functions obeying F* = F n *, where '*' denotes the complex conjugation. Then it is easy to show 
that eq. (17) yields, for the contour characterized by the abscissa c and the angle </> in fig. 2: 

1 f°° r t i 

F(x) = -J o dzlm [ eX p(i</>)x- c - zexp ^F n=c+zaLv{itl>) ] . (18) 

It is obvious from the discussion given above that the integral does not depend on c and <f>. 
However, for an efficient numerical evaluation a suitable choice of these parameters is very useful. 
For example, it is advantageous to choose </> > ir/2 in case F n is a known analytical function, 
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especially if this function does not fall off very rapidly for |n| — > oo. The factor exp (zlog ^ cos 
then introduces an exponential dampening of the integrand (which rapidly oscillates at small x) 
with increasing z, thereby allowing for a smaller upper limit z max in the numerical implementation 
of eq. (18). This procedure has been employed for the inversion of moments of parton densities 
and structure functions for the proton and the photon, e.g. in [6, 7] and [8], respectively. 

In general, however, the moments of the partonic cross section can only be calculated numer- 
ically using eq. (4), because no analytic continuation to small Ren, where the integral does not 
exist, is at our disposal. Likewise, in our case these moments are given by eq. (12) and do not 
behave uniformly for |n| — > oo. Especially, they grow exponentially along C\. Therefore, we will 
use the 'textbook contour' Co in the following and, with </> = 7r/2, eq. (18) simplifies to 

F(x) = - rdzRe\x- c - iz F n=c+iz ] . (19) 

TT JO 1 J 

We have applications in mind where F n = fi/ Pin h an , see eq. (15), and the numerical evaluation 
of the moments h an in eq. (12) is very time-consuming. Taking a different upper limit z max of 
the numerical z-integration or number of points for the integral at each step in the integration 
process is practically unfeasible in such a case. Instead, we want to fix z max at a value as small 
as possible in order to allow for an evaluation of eq. (19) with a rather small number of fixed 
moments. 

In this context, it is useful to consider that part of the ^-dependent inversion integrand, of 
which the analytical continuations are known together with the behaviour at large n, namely the 
parton densities and their evolution. Inspection of eq. (12) suggests that h an does not rise strongly 
with z along Co- Hence the large-n behaviour of the parton densities fi/ p>n can be employed to 
estimate the convergence of the complete integral in eq. (19) with F n given by eq. (15). A typical 
ansatz for the parton distribution functions of the proton at some reference scale fi^, denoted by 
/i/p(£)> is gi ven b Y [7, 9] 

e/i/p(o = ^ a (i-o*(i+7v£+--o • ( 20 ) 

The coefficients (3 can be estimated roughly by their counting rule values, e.g. (3 Y& \ 3, (3 g \ ue ~ 5. 
The Mellin transform of eq. (20) reads simply 

f i/Pin = A [B(a + n - 1,(3 + 1) + 1 B{a + n- 1/2,(3 + 1) + . . .] (21) 

with the Euler Beta-function B. If (3 is a positive integer, then this equation simplifies to 

fi/ P ,n = r — Uf , , ^ R u +■■■ = 0(l/n? +1 ) for n - oo . (22) 

[a + n — l){ct + n) . . . [a + n + (3 — 1) 

The evolution of these input moments is known analytically for arbitrary complex n [10], and 
the kernel Kij iTl (fi 2 , fi^) in 

fi/ P ,n(^ 2 ) = KijAv 2 1 Vo) fj/pAPo) ( 23 ) 
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generally leads to a slightly faster decrease of fi/ p>n for large n at fi 2 > fig. Hence a fall-off like 
1/n 4 can be safely used to estimate the practically required upper limit in eq. (19). 

For this purpose, we have numerically determined the upper limit z max sufficient to reach a 
1% accuracy of the Mellin inversion for the toy function F(x) = x~ x {l — a;) 3 . The results are 
displayed in table 1 for selected values of x. The rightmost pole is at n max = 1, and we have 
chosen c = 1.5. The larger c is, the more F n=c+ i z is flattened. Hence too large a value of c leads 
to an undesired rise of z max at small a;, where F n is integrated after multiplication with a rapidly 
oscillating function in eq. (19). Practically, c — n max ~ 0.5 works well for all ^-values of interest 
here, implying c ~ 1.8 for realistic small-£ parton densities [7, 9]. 



X 


10" 4 


10" 3 


0.01 


0.03 


0.1 


0.3 


0.8 


Zmax 


6.0 


5.5 


4.5 


4.0 


3.5 


5.0 


8.0 



Table 1: Upper limits z max numerically sufficient for a 1% accuracy of the Mellin inversion from 
eq. (19) for the function F(x) = x~ x {l — a;) 3 with c = 1.5. 

The integral eq. (19), truncated at z max , can now be performed by using a sufficiently large 
number of fixed support points, e.g. by a sum of 8-point Gaussian quadratures, see [11] for 
the weights and support points. In this way, everything except for the input-5-functions varied 
in a fit of the parton distribution functions is fixed; especially the time-consuming part of the 
kernels in eq. (23) and the moments h an of the partonic cross sections from eq. (12) can be 
determined once and then used unchanged in the calculation of physical observables for various 
parton densities. 

4 Application to Jet Physics at HERA 

To illustrate how the Mellin transform method can be used to fit the gluon density the gluon- 
induced (2+l)-jet cross sections were calculated in several bins for HERA energies of 820 GeV 
protons and 27.6 GeV electrons. Quark contributions were set to zero explicitly in the parton 
distribution function to reduce the number of moments needed for this case study. 

The program PROJET [12] based on the NLO matrix elements from [13, 14] was used, this 
allows to calculate jet cross sections in LO and NLO in the modified JADE scheme defined in 
the following way [15]: 

• Define a precluster of longitudinal momentum p r given by the missing longitudinal momen- 
tum of the event. 
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• Apply the JADE cluster algorithm [16] to the set of momenta {pi, . . . ,p n ,p r }, where 
pi,...,p n are the momenta of the hadrons visible in the detector. The resolution crite- 
rion is sa = 2piPj > VcutM 2 . Here M 2 is a mass scale and y cut is the resolution parameter. 

In the case of a theoretical calculation, p r is directly given by the momentum fraction of the 
proton not carried by the incident parton, and pi, . . . , p n are the momenta of the partons in the 
final state. In the following, we choose W 2 , the squared total hadronic energy, as the mass scale 
M 2 , since the proton remnant is included in the jet definition. 

The integration routine used in PRO JET is VEGAS [17, 18]. As is desirable for an exper- 
imental measurement, the phase space was binned in Q 2 and xb according to eq. (8); the bins 
are given in tables 2 and 3. In addition, the following typical HI detector cuts were applied, for 
which the motivation is explained in detail in [19]: 

• The invariant mass squared of the hadronic system W 2 was required to be larger than 
5000 GeV 2 . 

• The jet resolution cut y cut was set to 0.02. Lowering this value significantly below 0.01 
causes NLO corrections to dominate and leads to unphysical cross sections. It is important 
to note that £ > y cut as a consequence of the applied modified JADE algorithm. The region 
£ > 0.01 is however very interesting [20] for a precise determination of see also [21]. 

• The jets were required to lie in the polar angle range of 10° < 8j et < 145°. 

• For bins with Q 2 < 100 GeV 2 , the scattered electron had to have an energy of E\i > 14 GeV 
and the polar angle had to lie within the range of 160° < Q\i < 172.5°. 

• In the bins with Q 2 > 100 GeV 2 , the scaled photon energy y in the proton rest system had 
to be y < 0.7 and the scattered electron was required to have 10° < Q\i < 148°. 

In this list, angles and energies are defined in the laboratory frame, and angles are given with 
respect to the direction of the incoming proton. For each bin, 32 complex Mellin moments were 
calculated according to the prescription described in Section 2, cf. eq. (12). In all calculations, 
a s was computed to second order, and the NLO gluon distribution function of [7] was employed. 

A good convergence of the numerical calculations was found for c = 1.8, <f> = ir/2 and 
Zmax = 9, with a higher density of support points at lower z, as the influence is greatest there. 
For comparison, the cross section was also calculated directly, see eqs. (6), (7), (8). After inverting 
the product of the hard subprocess and evolved gluon density moments at the average Q 2 , the 
results were found to coincide at the per cent level. The detailed results can be found in tables 2 
and 3. In most bins, convergence was reached at z max = 3 (corresponding to 16 moments), the 
additional moments were used for safety. The convergence of the LO cross section was much faster 
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x B 


Q 2 [GeV 2 ] 


10" 4 


...1 


10" 3 


...1 


10" 2 


...1 


10" 1 


...1 


10 .. 


14 


62.80 


61.64 


28.09 


28.74 


— 


— 


— 


— 


14 .. 


18 


70.64 


69.72 


50.95 


49.97 


— 


— 


— 


— 


18 .. 


25 


85.82 


84.89 


71.03 


69.80 


— 


— 


— 


— 


25 .. 


40 


109.9 


108.8 


101.9 


101.1 


— 


— 


— 


— 


40 .. 


100 


— 


— 


123.8 


124.4 


14.51 


14.43 


— 


— 


100 .. 


300 


— 


— 


31.96 


32.18 


14.69 


14.76 


— 


— 


300 .. 


700 






28.97 


29.23 


25.18 


25.42 






700 .. 


. 4000 










10.22 


10.12 


0.96 


0.93 



Table 2: Comparison of cross sections with LO matrix elements 3 (in [pb]) obtained by integrating 
directly (left columns) or using the Mellin transform method (right columns). 





x B 


Q 2 [GeV 2 ] 


10" 4 


...1 


10" 3 


...1 


10" 2 


...1 


10" 1 


...1 


10 .. 


14 


58.48 


57.25 


26.60 


26.00 










14 .. 


18 


66.57 


65.90 


47.22 


46.69 










18 .. 


25 


82.48 


81.65 


67.99 


66.87 










25 .. 


40 


108.1 


107.4 


100.4 


99.71 










40 .. 


100 






126.1 


125.6 


14.07 


13.96 






100 .. 


300 






34.86 


34.52 


15.51 


15.31 






300 .. 


700 






31.34 


31.51 


27.01 


27.19 






700 .. 


. 4000 










11.18 


11.19 


0.99 


0.97 



Table 3: Comparison of NLO cross sections (in [pb]) obtained by integrating directly (left columns) 
or using the Mellin transform method (right columns). 

than in the NLO CclSCj clS for a given number of support points in VEGAS, the LO integration 
is more accurate due to the simpler integration kernel. The method works well for both LO and 
NLO. 

The number of points in the Monte Carlo integration was chosen such that the error returned 
by VEGAS was less than 1%. This number is, however, only a rough estimate [17, 18], and 
the achieved accuracy was studied by repeating the calculation for different random number 
generator seeds. The direct integrations performed here had a statistical variation of 2-3%. The 
partonic cross section from the Mellin transform method is implicitly integrated repeatedly by 
the calculation of the moments, which smoothes out statistical variations. The results were found 

3 Here, 'LO' means that the matrix elements were calculated in LO, but a s and the parton distribution functions 
in NLO to facilitate a comparison with the results of table 3. For a physically meaningful comparison of the LO 
with the NLO, a s and the parton distribution functions should be calculated in LO, if they are used in conjunction 
with the LO matrix elements. 
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to be more stable than the direct integration, which varied around the result obtained by the 
moment inversion. Even drastic errors of single moments or setting single moments to zero could 
be tolerated and led to a reproducible result. We conclude that this method is numerically very 
stable and that the accuracy is of the order of 1%. Increasing the accuracy requires increasing the 
number of support points for the integration, which would result in a dramatic increase in CPU 
time 4 . One has to keep in mind that an additional error source arises from the Mellin transform 
method, as for each experimental bin in x one has to calculate the difference of the cross sections 
depending on the bin boundaries in eq. (9), leading to error propagation. A strategy for bin 
optimization is under study. 



5 Summary 

We have outlined a new method, based on the Mellin transform for the fast evaluation of the 
convolution of a parton-level cross section with a parton density, which works for cross sections 
not showing a simple factorization behaviour. The method can be the basis for a fit of the gluon 
density from experimental data for the (2+l)-jet cross section in deeply inelastic electron-proton 
scattering, but it is also suitable for more complicated observables. It has been explicitly shown 
by a numerical study using realistic experimental cuts, that the method works in practice with 
a sufficiently high accuracy. 
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